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ABSTRACT 

The accretion-powered (non- X-ray burst) pulsations of XTE J1814-338 are modeled to determine 
neutron star parameters and their uncertainties. The model is a rotating circular hot spot and 
includes: (1) an isotropic blackbody spectral component; (2) an anisotropic Comptonized spectral 
component; (3) relativistic time-delays and light-bending; and (4) the oblate shape of the star due to 
rotation. This model is the simplest possible model that is consistent with the data. The resulting 
best-fit parameters of the model favor stiff equations of state, as can be seen from the 3-ct allowed 
regions in the mass-radius diagram. We analyzed all data combined from a 23 day period of the 2003 
outburst, and separately analyzed data from 2 days of the outburst. The allowed mass-radius region s 



for both cases only allow equations of state (EOS) that are stiffer than EOS APR (Akmal ct al. 



consistent with the large mass that has been inferred for the pulsar NGC 6440B (jFreire et al 



1998) 



2008) 



The stiff EOS inferred by this analysis is not compatible with the soft EOS inferred from a similar 
analysis of SAX J1808. 

Subject headings: stars: neutron — stars: rotation — X-rays: binaries — relativity — pulsars: indi- 
vidual: XTE J1814-338 



1. INTRODUCTION 

The accretion-powered millisecond-period X-ray pul- 
sars are promising targets for constraining the neu- 
tron star equation of state (EOS) through the mod- 
eling of emission from hot spots on the pulsar's 
surface. The fir st pulsar discovered in this class, 
SAX J1808.4-3658 (|Wiinands_fc van der Kliill998[ ). has 
a spectrum consistent (jGierlinski et al.ll2002h with emis- 
sion from a hot spot on the star's surface. Pulse shape 
modeling of rapidly rotating neutron stars relies on two 
relativistic effects: the gravitational bending of light rays 
reduces the modulation of the pulsed emission and de- 
pends on the mass to radius ratio Mj R; and the Doppler 
boosting due to the star's rotation creates an asymmetry 
in the pulse shape and depends on the star's radius R. 
These features, combined with reasonable models of the 
emission properties at the neutron star's surface can be 
used to constrain the neutron star's mass and radius and 
hence the EOS of supra-nuclear density matter. 

XTE J1814-338 (hereafter XT E J1814) was discov- 
ered during outburst in June 2003 (jMarkwardt fc Swank I 
2003), and is an accretion powered millisecond pulsar 
wit h spin frequency 314.36 Hz and orbital period of 4.3 
hr (jMarkwardt et al.ll2003l ). A detailed timing analy sis 
for XTE J1814 was performed bv lPapitto et al.l (|2007h to 
obtain accurate values for orbital period, projected semi- 
major axis, pulse spin frequency and spin down rate. A 
si milar analysis of the pulse ar rival times was car ried out 
bv lWatts fc Strohmaverl (|2006[ ) and[Cinmg| (|2007f) . which 
both included an analysis of phase lags. Soft lags were 
found in the 2-10 keV energy band, similar to those for 
SAX J1808-3658 and consistent with an origin in Doppler 
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boosting of a Comptonized pulse component with a much 
b roader emission pa tt ern th an the blackbody component. 

IStrohmaver et alj (|2003l ) found the same frequency 
in the X-ray bursts as was found in the persistent 
emission, but with a lower second harmonic content. 
I Watts et all ((2008) showed that the X-ray burst oscil- 
lations a re tightly phase-locked with the non-burst pul- 
sations. iBhattacharyya et al.l ((2005) modeled the oscil- 
lations during X-ray bursts with a hot spot model for a 
spherical star and for 2 equations of state. Using a large 
grid of models they found an upper limit on compactness 
Rs/R < 0.48, with Rs, the Schwarzschild radius. 

There are pulse shape models for a few other X-ray pul- 
sars. T he 1.2 s perio d X-ray pulsar Her X-l was mod- 
eled by ILeahv! (12001 ) using a model that includes ac- 
cretion columns. The model for Her X-l constra ins the 
neutr on star EOS to a fa ir ly mo derate stiffness (ILeahv 
2004). IZavlin fc Pavlovi (|1998t ) and iBogdanov eFal 
(2007) have modeled the X-ray emission from the 5.8 
ms period radio pulsar PSR J0437-4715 using a Hydro- 
gen atmosphere model . In the case of PSR J0437-4715, 
IBogdanov et alJ ((2007) found that a simple isotropic 
blackbo dy model is inco n sisten t with the data. In their 
models, IBogdanov et al.l (|2007l ) showed that the radius 
of PSR J0437-4715 must be larger than 6.7 km if the 
mass is 1.4M . Unfo rtunately the mass of t his pulsar is 
not well-constrained. iBogdanov et al.l (|2008[ ) have shown 
that constraints on radius for a number of other ms radio 
pulsars are also possible. 

Constraints on SAX J1 808. 4-3658 (with a spin period 
of 2.5 ms) were made by iPoutanen fc Gierliriskil (|2003l ) 
using da ta from the 1998 outburst. The modeling 
done by IPoutanen fc G icrlihski (|2003f ) included black- 
body emission from a hot spot that is Compton scattered 
by electrons above the hot spot. Their model makes use 
of a spherical model for the star's surface and does not 
include the effects of relative time-delays caused by the 
different time of flights for photons emitted from different 
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parts of th e star's surface . Mor e recently ICadeau et al.1 
(|2005[ ) and ICadeau et~afl (|2007| ) have shown that time- 
delays and the star's oblate shape are important factors 
that can affect the outcome of pulse-shape modeling for 
rapidly rotating pulsars such as SAX J1808. 4-3658 . The 
1998 outburst data for SAX J1808.4-3658 was revisited 
using a models that in cluded time-delays and oblate- 
ness (jLeahv et al.l I2008D with the result that the EOS 
for SAX J1808. 4-3658 is constrained to be very soft. 

In this paper we model the accretion-powered pulsa- 
tions of XTE J1814 using a hot-spot model. The hot- 
spot model allows for one or two circular hot spots with 
a two-component spectrum. The spectral model in- 
cludes isotropic blackbody emission and an anisotropic 
Compton-scattered component described by a power- 
law. The photons are propagated to the observer using 
the o blate Schwarzschild approximation (Mo rsink et al.l 
l2007f ) which allows the photon initial conditions to be 
placed on an oblate-shaped initial surface determined by 
an empirical formula. The Schwarzschild metric is used 
to compute the photon b ending angles and ti me delays 
since it has been shown (|Cadeau et al.l l2007t ) that the 
corrections induced by the Kerr black hole metric or a 
numerical metric for a rotating star are insignificant com- 
pared to the corrections induced by the oblate shape. In 
order to do the pulse-shape modeling, we construct light 
curves in two narrow energy bands, 2-3 keV and 7-9 keV. 
We first analyse a composite pulse-shape constructed 
from 23 days of data and then consider pulse-shapes con- 
structed from single days of data in order to determine 
whether variations of the pulse-shape with time are sig- 
nificant. 

The outline of this paper is as follows. In Section [2] the 
method used to construct the light curves and analyse 
them is outlined. The results of the best-fit models are 
presented in Section [3l A discussion of the results is 
presented in Section |H 

2. METHOD 

2.1. Construction of Light Curves 

Pulse shapes for the accretion-powered pulsations ar e 
constructed using the ephemeris given bv lChund((200l . 
Data is limited to the first 23 days of the 2003 outburst, 
June 5 to June 27 2003 (MJD 2452795-2452817), in order 
to avoid the later period of the outburst when the flux 
and pulse shape b ecame more erratic (see for example, 
iWatts et al.1 (|200l ). X -ray bursts were cut out of the 
data during the interval between 100 s before and after 
the start of each burst. 

Although the RXTE observations include data in the 
range of 2 - 50 keV, we have chosen to concentrate on the 
lower energy range from 2-10 keV for two reasons. First, 
the data is noisier at ener gies above 10 keV. S econd, the 
Chandra observations by iKrauss et al.l ((2005) constrain 
the spectrum in the 2-10 keV range. It is also useful 
to separate the data into narrow energy bands in order 
to separate the different spectral components. We have 
chosen two narrow bands, the 2-3 keV band and the 7-9 
keV band based on the Chandra spectrum. The narrow- 
band pulse shapes constructed using data from the full 
observation period (June 5 - 27) are shown in Figure [T] 

We have also investigated the variability of the pulse 
shape with time. In order to do this, the data was sep- 
arated into one-day segments and separate light curves 
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FlG. 1. — Pulse profiles for XTE J1814 constructed with data 
from all days between June 5-27 (excluding X-ray bursts). Light 
curves for two energy bands, 2-3 keV and 7-9 keV are shown. 
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FlG. 2. — One-day pulse profiles for XTE J1814 constructed with 
data from days June 20 and June 27 (MJD 2452810 and 2452817). 
Two energy bands are shown for each day. 

constructed for each day. It is not computationally feasi- 
ble to model all days simultaneously, so we instead focus 
on two separate days. The days were chosen by compar- 
ing light curves in the 2-10 keV range for different days 
using a x 2 test and selecting two days which differ the 
most from each other. This also has the effect of select- 
ing days with intrinsically smaller error bars. The days 
resulting from this selection process correspond to June 
20 and 27. Light curves for these two days, in the two 
narrow energy bands are shown in Figure O 

2.2. Analysis of Light Curves 

IKrauss et all (|2005l ) observed XTE J1814 with Chan- 
dra on June 20, 2003. They modeled the spectrum in 
the 0.5 - 10 keV range and found that the best fit solu- 
tion corresponds to a combination of a 0.95 keV black- 
body and powerlaw emission with a photon spectral in- 
dex of r = 1.4. The ratio of flux from the blackbody 
to the powerlaw in the ir model is about 10%. We use 
the IKrauss et al.l (|2005l ) spectral model and assume that 
it holds for the other days covered by the RXTE data. 
This assumption is motivated by the fact that the relative 
normalization of different energy bands is approximately 
constant from day to day, although the overall flux at all 
wavelengths changes wit h time. 

The spectral model of IKrauss et al.1 (|2005h motivates 
the use of two narrow bands in our pulse shape models. 
A low energy band is necessary in order to capture the 
blackbody component, so we choose the lowest possible 
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XTE energy band at 2 - 3 keV. The spectrum in this 
band is dominated by the powerlaw component, but the 
blackbody contribution is still important. We also choose 
the 7-9 keV band as the highest energy band covered 
by the Chandra observation. In this high energy band 
the blackbody flux is negligible. 

Our method for modeling the obse rved emission is ver y 
similar to the method presented in iLeahv et al.1 $2008). 
The spectral model has three components: (1) Comp- 
tonized flux in the high energy band (7-9 keV); (2) Comp- 
tonized flux in the low energy band (2-3 keV); and (3) 
blackbody flux in the low energy band. The observed 
flux for the ith component, Fi, integrated over the ap- 
propriate observed energy band is given by 

Fi(E) =/ 4 77 3+r *(l-a^). (1) 

In equation ([!]). Jj is a constant amplitude, r\ is the 
Doppler boost factor, I\ is the photon spectral index in 
the star's rest frame, [i is the cosine of the angle between 
the normal to the star's surface and the initial photon di- 
rection, and the constant describes the anisotropy of 
the emitted light. For a definition of r\ as well as a more 
complete descriptio n of the modeling method, please see 
ILeahv et al l $20081) . 

In our modeling, the amplitudes I\ and I2 are free pa- 
rameters while the_third amplitude I3 is defined through 
the constant b = F 3 /F 2 , the ratio of the phase-averaged 
blackbody to Comptonize d flux in the low-ene rgy band. 
In the spectral model bv iKrauss et al.1 $2005), b = 0.1, 
but we include this parameter as a fitting parameter with 
lcr limits from their spectral model. The photon spec- 
tral indices for the Comptonized components are fixed at 
Ti = T2 = 1.4 as given by their model. In the narrow 
range of the low-energy band the blackbody component 
of 0.95 keV can be modeled by a powerlaw with photon 
spectral index = 0.85. The anisotropy parameters 
for the Comptonized components a\ = a 2 = a are as- 
sumed to be equal, and the parameter a is kept as a free 
parameter. 

In the modelling of the non-accreting ms pulsars, it 
was f ound (|Zavlin fe Pavlov|[T998t iBogdanov et al.ll2007l 
2008) that a limb-darkened Hydrogen atmosphere spec- 
tral model is required by the data. It is reasonable to ex- 
pect that that the blackbody component of the spectrum 
should also be limb-darkened. We tested this hypothesis 
by multiplying the blackbody flux by a limb-darkening 
function of the form e _T / p . We then computed the bestfit 
neutron star models for two type of models: (1) models 
with non-zero optical depth r and (2) models with zero 
optical depth. The bestfit models for these two cases 
are almost identical: the mass and radius of the bestfit 
model changes by less than 0.5% when a nonzero opti- 
cal depth is added, and the value of 6\ 2 =0.1 when 
the limb-darkening is added. Since the change in \ 2 and 
the physical parameters are negligible we conclude that 
adding an extra parameter to model limb-darkening is 
not warranted by the data. The reason for this is due to 
the Chandra model which restricts the blackbody con- 
tribution in the 2-3 keV band to only 10% of the Comp- 
tonized contribution, and effectivly sets the blackbody 
component to zero in the high energy band. Since the 
Comptonized flux is dominant and has fan-beaming in- 
cluded, small changes to the anisotropy of the blackbody 
component don't affect the final models. For this reason 



we have set the anisotropy parameter for the blackbody 
component to zero (0,3 = 0). This is consistent with the 
results found for SAX J1808.4-3658 $Leahv et al.ll200l 
which also did not require any limb-darkening. The final 
set of free parameters describing the spectrum are I\ , 1% , 
b and a. 

In order to fit a set of light curves we also need to 
introduce a set of parameters describing the star and the 
emission geometry. These parameters are the mass M 
and equatorial radius R of the star, the co-latitude of 
the spot 9, the inclination angle i as well as a free phase 
4>. The radius of the circular spot (in the star's rest 
frame) is kept fixed at 1.5 km, as g iven by the Chandra 
spectral model (|Krauss et al.l 120051 ). 

Our models make use of light curves for two different 
days' data, which requires a separate set of parameters 
for each day. However, on the two different days the pa- 
rameters M , R and i do not change. In order to simplify 
the analysis, we also assume that the photon spectral 
indices and the parameters a and b are also fixed. The 
full set of free parameters are: {Ii, 1 2 , 0, 4>} for each day 
plus M, i?, i, a and 6, for a total of 13 free parame- 
ters. However, for each of our fits, the ratio M/R is kept 
fixed, so for any one value of M/R, there are only 12 free 
parameters. 

We use the obl ate Schwarzschild approximation 
(jMorsink et al.ll2007h to connect photons emitted at the 
star's surface with those detected by the observer. In pre- 
vious studies (jCadeau et alJl2005l [2007]) we have shown 
that, to the accuracy required for extracting the param- 
eters of a rapidly rotating neutron star, it is sufficient to 
use the Schwarzschild metric to compute the bending of 
light rays and the relative time delays of photons emitted 
at different locations on the star. The extra time delays 
and light bending caused by frame-dragging or higher 
order rotational corrections in the metric are negligible. 
However, the rotation of the star causes a deformation 
of the star into an oblate shape, which changes (relative 
to a sphere) the directions that photons can be emit- 
ted into. We have developed a simple approximation 
(jMorsink et al.ll2007| ) that allows an empirical fit to the 
oblate shape of a rotating star to be embedded in the 
Schwarzschild metric and make use of it in this analysis. 

3. RESULTS 

3.1. Evidence for a Second Spot 

The pulse profiles in Figure [T] show a feature in the 
phase interval between 0.24 and 0.4. This feature is seen 
in all of the other energy bands as well. In order to inves- 
tigate the nature of this feature, we restrict the analysis 
here to just the 7-9 keV light curve. Since the blackbody 
contribution in this energy band is negligible, we use a 
simplified model which only includes the Comptonized 
component of the radiation. 

The simplest model for the emission is a single spot. 
We fitted the 7-9 keV light curve shown in Figure[3]with a 
single spot model by first fixing 2M/R = 0.4 and varying 
the parameters M,i,9,a,I and </>. (Similar results are 
obtained for other values of 2M/R.) Since there are 32 
data points this corresponds to 25 degrees of freedom. 
The best fit solution for a single spot model has x 2 = 
50.7, which is not a very good fit. This best-fit solution 
is shown as a solid curve in Figure [3] 

We now turn to a two-spot model, where the second 
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Fig. 3. — Models used to fit the data in the 7-9 keV band (data 
combined from all 23 days). A one-spot model that uses all data 
points (both squares and circles) results in the best fit solid curve 
with x 2 /dof = 50.7/25. A two-spot model that uses all data points 
results in the best fit dotted curve with x 2 /dof = 28.3/23. A one- 
spot model that omits the circle bins results in the best fit dashed 
curve with x 2 /dof = 19.2/19. 

spot is allowed to have an arbitrary location relative to 
the first spot, but the spectrum is assumed to be the 
same. The introduction of a second spot introduces three 
new parameters to the model: an intensity and two an- 
gles. The resulting fit for 2M/R = 0.4 has \ 2 = 28.3 
for 23 degrees of freedom, which is a significant improve- 
ment. The light curve for this model is shown with a 
dotted curve in Figure [3] The mass for this model is 
2.04 Mq. In this model, the second spot's location is sit- 
uated so that the second spot is almost never seen and its 
light only contributes during the phase interval between 
0.24 and 0.4. Outside of this interval the light received 
only originates from the primary spot. This suggests a 
simpler one-spot model where bins in the phase interval 
0.24 to 0.4 (marked with a circle in Figure[3]) are removed 
from the data set. This reduces the degrees of freedom 
to 19 (32-6 data points and 5 parameters for a one-spot 
model). The resulting best-fit model for 2M/R = 0.4 
has x 2 = 19.2, which is also a significant improvement 
from the one-spot model that uses all of the data. The 
mass for this model is 2.08M Q and the light curve for 
this model is shown as a dashed line in Figure [3l 

Comparing the two-spot model with the one-spot 
model with the circle-bins excluded, we see that the dif- 
ference in \ 2 is not significant, and there is little change 
between the best-fit values of mass and radius. This leads 
us to the conclusion that there is good evidence for a sec- 
ond spot (or a feature that mimics a second spot), but 
that the amount of data encoded in those bins affected 
by the second spot is not sufficient to allow us to model 
the details of the second spot with any confidence. Since 
the inclusion or exclusion of the second spot does not 
change the best-fit values of the main physical parame- 
ters of the star (M, R, i) it is more appropriate to choose 
the simpler one-spot model. The qualitative results do 
not change when we look at different energy bands. As a 
result, for the remaining modeling reported in this paper 
we only use one-spot models where data in the 0.24 - 0.4 
phase range is removed from the analysis. 

3.2. Best-fit Models using Data from All Days 

Our procedure for modeling the two-energy band data 
shown in Figure [T] is to assume a one-spot model that in- 
cludes both blackbody and Comptonized emission. Data 



TABLE 1 
Comparison of the Two Minima. 



Model 


2M/R 


M 


R 


e 


i 


a 


M c 


x 2 /dof 








km 


deg. 


deg. 




Mq 




High Mass 


0.3 


2.86 


28.7 


66.8 


11.8 


0.81 


1.70 


55.9/42 


Low Mass 


0.3 


1.95 


19.8 


42.4 


25.0 


0.59 


0.55 


61.0/42 



in the phase period between 0.24 and 0.4 is omitted, as 
described in section 13.11 We do a number of fits, each 
with a fixed value of 2M /R. (Fixing the ratio of M/R 
simplifies the fitting procedure since the light-bending 
and time delays depend on this ratio.) Once 2M/R is 
fixed, all other parameters are allowed to vary and the 
minimum value of \ 2 is found. In addition to the pa- 
rameters described in Section |2~2"1 we also added two pa- 
rameters corresponding to DC offsets for the two energy 
bands. This allows for small errors in the background 
subtraction. Once 2M/R is fixed, we have a total of 
10 free parameters: M, 8, i, a, b, 2 amplitudes, 2 DC 
offsets and one overall phase. (The parameter b is re- 
stricted t o have a value that is within 1 a of the value 
found bv iKrauss et alJ 1)20051) .) Since each energy band 
has 32 points, but we exclude 6 of these points we have 
64— 12 — 10 = 42 degrees of freedom. 

For a fixed value of 2M/R, we find that there are two 
local minima. These two minima are shown in Table [T] 
and we label these two best-fit solutions as the high and 
low mass solutions. The lowest value of x 2 corresponds 
to the high mass solution and the lower mass solution has 
a higher value of x 2 ■ Although the high mass solution is 
a better fit, we exclude this solution on physical grounds. 
First of all, it requires a neutron star radius of 28 km, 
which is not allowed by any known equation of state. 
Secondly, once the neutron star mass and the inclination 
angle are known, the companion's mass can be calculated 
(shown in the column labeled M c in Table [T]) . In the 
high mass case the companion's m ass is 1.7 M Pl . Due to 
the dim nature of the companion, IKrauss et ail (|2005D 
have shown that the companion (if a main sequence star) 
would have to have a mass that is no bigger than 0.5M Q . 
Clearly this excludes the high mass solution but allows 
the lower mass solution. For these reasons, we exclude 
the high mass solutions. 

In Table [2] the best-fit solution for each value of 2M/R 
is shown. In each case only the lower mass solution is 
shown. The best-fit solution shown as a solid curve in 
Figure [1] corresponds to the 2M/R = 0.4 solution shown 
in the Table [2j Although we call this these solutions 
"lower mass" , clearly they still correspond to high mass 
neutron stars. In Table[2j only solutions for 2M / R = 0.3 
to 0.6 are shown. In the case of less compact stars, such 
as 2M/R = 0.2, the only solution corresponded to the 
"high mass" branch of unphysical solutions. We did not 
test solutions that were more compact than 2M /R = 0.6 
since these solutions would allow spots to have multiple 
images, and our program is unable to handle multiple im- 
ages. Technically, the solutions with 2M/R = 0.6 could 
have spots with multiple images, but we have checked 
that the relative values of 9 and i do not lead to this 
problem for our solutions for the most compact stars. 

In Figure 2] we show a number of mass-radius curves 
for stars spinning at 314 Hz as well as the 2- and 3- 
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TABLE 2 

Best-fit solutions using data from all days. 



2M/R 


M 


R 


e 


i 


a 


M c 


X 2 /dof 




M 


km 


deg. 


deg. 




Mq 




0.3 


1.95 


19.8 


42.4 


25.0 


0.59 


0.55 


61.0/42 


0.4 


2.45 


18.4 


47.0 


24.2 


0.61 


0.65 


61.2/42 


0.5 


2.38 


14.2 


36.7 


39.4 


0.59 


0.39 


61.9/42 


0.6 


2.42 


11.9 


40.9 


42.8 


0.59 


0.37 


62.6/42 
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Fig. 4. — Best-fit mass and radius values using the combined data 
from all days, separated into two narrow energy bands. Contours 
shown are for 2- and 3-<r confidence levels. Mass-Radius curves for 
stars spinn ing at 314 Hz are sh own as s olid curves. The EOS shown 
ar e: APR fAkmal et al.lll99Sti, BBB1 lIBaldo et a l. 1997), ABPR2- 
3 l|Alford et al.H2005l ). H3-7 ILackev Nawar. fc OwerJH2006fl an d L 
(mean- field theory, pure neutrons l|Pandh*arl^mde~et£d!lll976l \). 

a confidence regions for the "low mass" solutions. (In 
this figure, the radius R refers to the equatorial radius.) 
These regions are found by fixing the value of 2M/R, 
and varying all other parameters and finding the solu- 
tions that have \ 2 larger than the global minimum of 
x 2 nm = 61.0 by Sx 2 = 4 (2 a) or 9 (3 a). The region 
allowed with 99.7% confidence (3 a) only includes large 
stars with high mass. Of the equations of state displayed 
in FigureUl the only one lying in the 3 a allowed region is 
L, corresponding to pure n eutron matter computed in a 
mean field approximation (|Pandharipande et alj Il976f ) . 
A pure neutron core is unlikely to be the correct de- 
scription of supra-nuclear density matter. However it is 
possible for an EOS that includes some softening due to 
the presence of other species to be allowed by this data. 

3.3. Best-fit Models Using Data From 2 Days 

It is possible for variability of the data to affect the 
fit results, as app ears to be the case for SAX J1808 
(|Leahv et alJ 12008). For this reason we have rebinned 
the data into one-day segments in order to see if there is 
any significant change in the pulse shape on a day to day 
basis. We performed a x 2 test to see how closely each 
day's data matched the other days' data. Comparisons 
of one day with an adjacent day gave values of x 2 /dof 
ranging from 0.8 to 1.6, indicating day-to-day changes 
are small. The largest change is between day 810 (June 
20) and day 817 (June 27) with x 2 /dof = 4.8. Light 
curves in the two narrow energy bands for these two days 
are shown in Figure O 

The data corresponding to these two days is binned 
into 32 bins per period. Since we continue to remove the 
6 bins corresponding to the second "spot" , this corre- 



TABLE 3 

Best-fit solutions using data only from Days 810 &; 817. 



2M/R M R e 1 8 2 i a M c x 2 / dof 

Mq km deg. deg. deg. Mq 

0.2 1.18 18.4 33.6 34.7 33.8 0.55 0.29 123.9/92 

0.3 1.71 17.3 34.1 35.3 32.0 0.55 0.39 124.6/92 

0.4 2.13 16.1 35.1 36.6 33.6 0.55 0.43 125.4/92 

0.5 2.41 14.4 37.1 38.7 36.3 0.55 0.43 127.5/92 

0.6 2.55 12.6 39.5 41.4 39.9 0.56 0.41 132.0/92 
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Fig. 5. — Best- fit mass and radius values using data only from day 
810 and 817, separated into two narrow energy bands. Contours 
shown are for 2- and 3-cr confidence levels. EOS labels are the same 
as in Figure [4] 

sponds to a total of 104 data points. We fit the data 
for these two days by assuming that the parameters 
M,R,i,a are the same for both days. The spot's lat- 
itude is allowed to vary, as are the amplitudes of the 
energy bands and the value of phase. Since the DC con- 
tributions found in our previous fits are very small, we do 
not include terms for DC offsets. In additi on, we keep the 
value of b fixed at the iKrauss et al.l (|2005f ) value in order 
to simplify the fits. This corresponds to 12 parameters, 
and a total of 92 degrees of freedom. 

The best-fit solutions for this 2-day joint fit are shown 
in Table [3) For each fixed value of 2M/R we only found 
one minimum, unlike the case with all days included. 
The angular locations of the spot on the two days are 
are labeled Q\ and 6*2- The change in angular location of 
the spot between the two days is less than 2° in all cases. 
The solutions continue to have large masses and radii, 
as in the case of fits using all of the data. In the case 
of 2M/R = 0.2 a low mass (1.2 Mq) solution is allowed, 
but it has a very large radius. The 2- and 3-er confidence 
regions for the two-day joint fits arc shown in Figure O 
The 3-cr confidence region is somewhat larger than the 
same region computed using all of the data, but the two 
methods have a significant overlap. The 2-day joint fit 
also only allows the stiffest EOS L. 

3.4. Dependence of Models on Assumed Parameters 

The models in th is paper depe n d on t he results of the 
spectral models of IKrauss et alJ (j2005| ). We now con- 
sider the effect of allowing the parameters in the spec- 
tral model to vary within the error bars. As mentioned 
in Section 12.21 we already allow the ratio of the black- 
body to powerlaw c omponents (b ) to va ry within the la 
limits given by the IKrauss et al.l (|2005l ) spectral model. 
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Another spectral parameter that could affect the fits is 
the photon spectr al index T for t he pow er law component 
of the spectrum. iKrauss et alj (2005) found a value of 
r = 1.41 ± 0.06 and in all of our models presented in the 
previous section we kept the photon spectral index fixed 
at a value of T = 1.40. We would expect that a larger 
value of r would allow for smaller stars. This is because 
the flux in Equation {T]) is proportional to the Doppler 
boost factor r\ raised to the power T + 3. The Doppler 
boost factor is mainly responsible for introducing higher 
harmonics into the signal, so a larger value of T creates a 
more asymmetric pulse shape. In order to compensate, 
the best-fit solution will require a smaller value of stellar 
radius R in order to decrease the value of r\. In order 
to test the dependence of the best-fit values of the pa- 
rameters on T, we chose a value of T = 1.50 which is 
somew hat larger than the range allowed bv lKrauss et al.1 
( 2005]) and fit the data using the same method described 
earlier in this paper. The results of the best-fit parame- 
ters for the two values of T for the case of 2M/R = 0.4 
are shown in Table 0J As expected, increasing T allowed 
for a smaller star, but the decrease is only by 3%. Sim- 
ilar results occur for other values of M/R. Clearly the 
dependence on the photon spectral index is not sensitive 
enough to affect the resulting large size of the best-fit 
stars. 

In our models we keep the spot size (as measured on the 
star's surf ace) fixed at a diam eter of 3 km. In previous 
modeling dLeahv et all l2008h of SAX J1808.4-3658 we 
found that the final values of the best-fit parameters 
were not sensitive to changes in the spot size, assum- 
ing that the spot is small compared to the size of the 
star. For this reason we have kept the spot size fixed 
at 3 km for all models in this analysis. We have also 
assumed the hot spot is circular, alth ough recent MHD 
models (jKulkarni fc Romanovall2005h more complicated 
spot shapes. In this paper we have only attempted to 
use the simplest possible model that is still consistent 
with the data. Adding extra parameters to our models 
in order to describe more complicated spot shapes is not 
yet warranted by the quality of the data. 

For all models computed so far, we have made use 
of an empirical formula for the oblate shape of the 
star, and have included the relative time-delays for pho- 
tons emitted from different parts of the star. In Table 
[4] the effects of oblateness and time-delays on the fits 
are shown. For the model labeled "sphere", a spheri- 
cal initial surface was assumed, but relative time-delays 
were included in the computation. The resulting best- 
fit solution is about 10% larger than the corresponding 
oblate model (labeled "oblate" and V = 1.4 in Table 
0J. This shrinkage of the star's radius when oblate- 
ness is included ha s been observed in the modeling of 
SAX ,11808.4-3658 dLeahv et all 12001 . For the model 
labeled "no td" time-delays were omitted from the calcu- 
lation and a spherical surface was used. Comparison of 
the two spherical models in Table[4]shows that the model 
that includes time-delays is about 3% smaller than the 
model that omits time-delays. 

3.5. Comparison with X-ray Burst Data 

In our analysis of XTE J1814 we have only included 
data from the accretion-powered pulsations and have 
omitted any data corresponding to an X-ray burst. 



TABLE 4 

Dependence of Models on Parameters. Joint fits for two energy 
bands for two separate days (810 and 817). 



Model 


r 


2M 
R 


M 

Mq 


R 
km 


01 

deg. 


e 2 

deg. 


i 

deg. 


a 


M c 

Mq 


x 2 /dof 


oblate 


1.4 


0.4 


2.13 


16.1 


35.1 


36.6 


33.6 


0.55 


0.43 


125.4/92 


oblate 


1.5 


0.4 


2.07 


15.6 


34.5 


35.9 


34.3 


0.55 


0.41 


126.6/92 


sphere 


1.4 


0.4 


2.38 


17.6 


32.8 


34.3 


31.6 


0.54 


0.49 


124.7/92 


no td 


1.4 


0.4 


2.45 


18.5 


49.0 


50.9 


21.1 


0.60 


0.76 


126.8/92 



iBhattacharvva et al.1 (|2005t ) analyzed the light curves 
constructed from the X-ray bursts for this neutron star. 
In their analysis they assumed a spherical surface for 
the star and traced the paths of the X-rays using the 
Kerr metric. They also made use of a limb-darkened 
blackbody emission (2 keV) spectral model appropri- 
ate for X-ray bursts. Due to the method that they 
adopted, it was necessary for them to assume one of two 
different equations of stat e. The stiffer EOS used by 
Bhattach arvva et al.l (|2005[ ) is the same as the EOS that 
we label "APR" an d is the A18+fo;+UrX computed by 
Akmal et al.1 (119981). The a nalysis of the X-ray bursts by 



Bhattacharvva et all (|2005l ) allows the APR EOS, while 



our analysis of the accretion-powered pulsations only al- 
lows stiffer EOS. From their analysis it is difficult to de- 
termine whether or not their analysis of the X-ray burst 
data is consistent with a very stiff EOS, as indicated by 
o ur analysis. 

IWatts et al.l (|2005f ) provide a detailed analysis of many 
aspects of both the X-ray bursts and the non-burst emis- 
sion. One of the quantities that they measured was 
the fractional amplitude of the pulsations at the funda- 
mental frequency and the first harmonic for both the 
burst and non-burst emission. They found that the 
non-burst emission (modeled in this paper) has a larger 
harmonic content than the burst emission studied by 
IBhattacharvva et al.1 (|2005h . Since Doppler boosting is 
partially responsible for the harmonic content, it is per- 
haps not surprising that our analysis of the non-burst 
pulsations implies a larger Doppler factor for the star, 
which in turn implies a larger radius than a similar anal- 
ysis for the X-ray burst oscillations. 

4. DISCUSSION 

We have analyzed the accretion-powered (non- X-ray 
burst) pulsations of XTE J1814 using a hot spot model. 
Our modeling includes (1) an isotropic blackbody spec- 
tral component; (2) an anisotropic Comptonized compo- 
nent; (3) relativistic time-delays; (4) the oblate shape of 
the star due to rotation. The model presented in this 
paper is the simplest possible model that is consistent 
with the data. The resulting best-fit models favor stiff 
equations of state, as can be seen from the 3-cr allowed 
regions in Figures 2] and [5] In Figure 2] all data from a 
23 day period of the 2003 outburst were included, while 
for Figure [5] data from only 2 days were included. The 
allowed regions for the two data sets differ slightly, but 
both only allow equations of state that are stiffer than 
EOS APR (| Akmal et al.lfl99l . 

It is interesting that a la rge mass has been inferred 
for the pulsar NGC 6440B (jFreire et al.ll2008f) through 
measurements of periastron precession. Assuming that 
the observed periastron precession is purely from rela- 
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tivistic effects, t he pulsar's mass is M = 2.74 ± 0.21M Q 
(l-o" error bars) (jFreire et alj|200g ). If the mass really is 
this high, it would be consistent with the stiff equations 
of state allowed by our analysis of XTE J1814 . How- 
ever, it is still possible that the large periastron preces- 
sion observed for NGC 644 0B could be caused by a very 
rapidly rotating companion iFreire et al.l (|2008f), in which 
case the pulsar's mass would be smaller and compatible 
with more moderate equations of state. A high mass for 
SAX J1808. 4-3658 is a lso inferred by obs ervations dur- 
ing its quiescent state dHeinke et al.ll2007T) . Mode lling of 
the neutron star X7 in 47 Tuc (jHeinke et alll2006h allows 
for a high mass neutron star, although for X7 a low mass 
neutron star is also allowed. 

Similar hot spot models of SAX J1808.4-3658 imply a 
soft equation of state and a c olumn model for H er X- 
1 also implie s a soft EOS (see iLeahv et al.l (|2008l ) and 
iLeahvl (|2004h for details). The best-fit pulse-shape mod- 
els found for XTE J1814 have mass and radius incom- 
patible with the 3-cr allowed regions of mass and radius 
found for SAX J1808.4-3658 or Her X-l. Some possible 
interpretations could be (1) time- variations in the pulse 
profile of SAX J1808. 4-3658 led to an underestimate of 
the star's radius; (2) the equation of state for dense mat- 
ter has a two-phase nature allowing both large and small 
compact stars; or (3) the simple hot-spot model doesn't 
describe one or more of these pulsars. 

The first reason is a factor for SAX J 1808.4-3658 . It 
was demonstrated in lLeahv et al.l ()2008l ) that the analy- 
sis of the pulse profile averaged over a long observation 
gave significantly different results than the pulse profile 
obtained from a much shorter observati on. Supporting 
this c onclusion is a recent analysis by lHartma n et al.l 
(2008) of data from the 1998, 2002 and 2005 outbursts 
of SAX J1808. 4-3658 showing a great de al of variation in 
the pu lse shape over time (see Figure 3 of Hartman et all 
(2008)). The pulse-shape analysis bv ILeahv etafl 1)20081 ) 
only made use of the data from the 1998 outburst. The 
1998 data is very sinusoidal in na ture and has very little 
harmonic content. The results of lHartman et all (2008) 
show that the later outbursts have a stronger harmonic 
content. Since a larger radius star can produce a stronger 
harmonic content, it is possible that the addition of 
data from 2002 an d 2005 will alter the conclusions of 
ILeahv et all (|2008l ) about SAX J 1808. But the effect of 
pulse shape variability is not a factor for Her X-l where 
the pulse shape has high stability. 
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The second reason above, i.e. a bimodal equation 
of state, is a possible, but speculative, solution to the 
greatly different allowed regions for M and R. In this sce- 
nario, there is still only one baryonic equation of state 
and one quark matter EOS, but above a certain criti- 
cal density, p cr it , the whole star makes a transition from 
baryonic matter to quark matter. Then for stars with 
central density p c having p c < p cr u, the M vs. R rela- 
tion follows a stiff baryonic EOS, somewhat like EOS L, 
whereas for p c > p cr it the star has converted to a quark 
star and lies on a quark matter EOS curve in the mass- 
radius diagram. In the case of EOS L, the 3<r region with 
M < 2.7Mq would require a value of p cr i t ~ 10 15 g/cm 3 . 
Since quark matter EOS curves have a lower maximum 
mass than baryonic EOS curves, any baryonic star that 
makes the transition and has mass above the quark star 
maximum mass must lose mass to end up as a stable 
quark star. The mass loss depends on the physics of the 
transition process and is likely vary from star to star. 
In this scenario, we interpret XTE J1814 to lie on the 
baryonic branch of the mass vs. radius diagram and 
SAX J1808. 4-3658 and Her X-l to lie on the quark mat- 
ter branch. 

The third listed reason for the discrepancy, that the 
emission region models are too simple to represent the 
actual emission regions on the stars, is a definite possi- 
bility. For instance, if the emission is coming from the 
magnetosphere, then our models are incorrect. Alterna- 
tively, the emission may arise from surface spots, but the 
region's shape might be more complicated than a circle. 
This can only be tested by constructing more complex 
models and applying them to the observed pulse shapes. 
However with a more complex model with more parame- 
ters describing the model, better pulse shape data is re- 
quired to constrain the model parameters. Future work 
is planned to explore more complex emission models, and 
to test whether these resolve the apparent discrepancy in 
mass and radius for different pulsars. 
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